Critical Velocity of Superfluid Flow past Large Obstacles in Bose-Condensates 
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By considering the stability of potential flow of a superfluid around large obstacles of size R, we 
derive an analytical result for the critical velocity which is of order v c ~ h/mR, scaling inversely with 
obstacle size, in contrast to what is obtained from a Landau criterion. Our results are compared with 
numerical solutions of the Gross-Pitaevskii equation and with recent measurements of the critical 
velocity in Bose- Einstein condensates of dilute atomic gases. 
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Dissipationless flow for small, finite velocities is a basic 
and, originally, defining property of superfluids JjJ. For 
Bose-Einstein-condensates in dilute atomic gases, this 
phenomenon was recently investigated by stirring the 
condensate with a blue detuned laser beam ||, ||. The 
associated time dependent repulsive potential of macro- 
scopic size R, much larger than the healing length £, gives 
rise to a finite dissipation only above a critical velocity 
v c , where superfluidity breaks down. The rate of dissipa- 
tion can be obtained either by measuring the amount of 
heating via the resulting depletion of the condensate 
or, more directly, by observing the asymmetry in den- 
sity associated with a finite pressure difference across 
the moving object B. The latter phenomenon is in fact 
completely analogous to the current induced asymmetry 
around defects in a Fermi liquid, as predicted long ago 
by Landaucr as the basic microscopic origin of residual 
resistance §, |. The measured values of v c are found 
to be a small fraction v c « (0.1 — 0.25)c s of the sound 
velocity at the trap center with a numerical factor, which 
appears to increase with density (|] . These results are in 
rough agreement with numerical solutions of the Gross- 
Pitaevskii equation (GPE) for a moving cylindrical per- 
turbation in a homogeneous B, fjj or inhomogeneous ||] 
condensate. In the homogeneous case, for a fixed cylinder 
size, the critical velocity is around 0.45c s . Qualitatively 
this result may be understood in the spirit of the Landau 
criterion [Q, by arguing that there is no dissipation as 
long as the local velocity u(x) at any point remains be- 
low c s . For an incompressible flow, u(x) has its maximal 
value, which is twice the object speed, at the side of the 
cylinder \x\ = R, 6 = ±n/2 ||. Since the local sound 
velocity is reduced by the depletion of the actual conden- 
sate density close to the obstacle, a value v c smaller than 
c s /2 (or 2/3 c s for a sphere) but independent of R is ex- 
pected from this type of argument. In an inhomogeneous 
situation like that for atomic condensates in a harmonic 
trap, the critical velocity will be lowered further because 
the density decreases away from the trap center. 

Our aim in the present work is to derive an analytical 
criterion for the critical velocity of superfluid flow around 
macroscopic obstacles of size R ^> £ by considering the 
linear stability of the dissipationless potential flow below 
v c towards the generation of vortices. For a strongly 
repulsive potential, v c is of order bh/mR i.e. independent 



of the density and inversely proportional to the obstacle 
size. Since h/m£ = \/2c s , the critical velocity coincides 
with a fraction of the sound velocity for obstacle sizes 
R w (10— 20)£, which are in fact typical values considered 
both experimentally and in the numerical simulations. 
As pointed out by Nozieres and Pines [jj], the prediction 
v c ~ h/mR is a generic result for the critical velocity 
due to vortex generation, as shown e.g. by Feynman's 
estimate v c w In § <C c s jl0| for superfluid flow in a 
long channel with diameter D ^> £ |fu]| . A very similar 
result for a cylindrical trap was obtained very recently 
by Crescimanno et al p2| . 

The dynamics of a dilute, weakly interacting Bose- 
Einstein-condensate (BEC) is accurately described by the 
GPE for the complex condensate wave function <f> |L3| |. 
Since we are considering a bulk situation, it is convenient 
to normalize $ according to J \&\ 2 d 3 x — N. In the rest 
frame of the obstacle, where the problem is stationary, 
the energy functional J-[$>] of a homogeneous condensate 
moving with velocity v, has the standard form 



Here 
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is the condensate energy functional in its rest frame and 

= J d d x^ [$V$* - $*V$] (3) 

is the momentum functional of the moving condensate. 
For dilute gases, the interaction parameter g = Anf^a/m 
is determined by the scattering length a, which is as- 
sumed to be positive. Splitting the condensate wave 
function Q(x) — f{x)e lip ^ into magnitude and phase, 
the energy functional ([[]) takes the form 
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In order to calculate the critical velocity, we determine 
the limit beyond which small fluctuations around the dis- 
sipationless potential flow below v c are no longer stable. 
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Since rottt vanishes identically below v c , the condensate 
velocity relative to the obstacle at rest can quite generally 
be written in the form 

u(x) = — Vip(x) + v , (5) 
m 

where the term li/raV^ accounts for the backflow in- 
duced by the presence of the obstacle. Far from the ob- 
stacle, where the density f 2 is constant, the backflow is 
necessarily of a dipolar form 



h , (dn(v ■ n) — v) 

-Vip = - V — ^ — ^ '- 
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since the dipole is the only function obeying both con- 
ditions of vanishing vorticity rotw = and stationarity 
= div/ 2 it w / 2 dfvu. Here n is the unit vector in the 
direction of x and d — 2,3 the relevant dimensionality 
(d = 2 for a moving cylinder, d = 3 for a sphere). For a 
weak obstacle potential Vb <C Mi the response of the su- 
pcrfluid to the movin g ob ject can be treated in first order 
perturbation theory [Q. The resulting dipole strength 
Ad ^ ~ • V (q = 0) is determined by the bulk compress- 
ibility Tin = i I 2 = (mc 2 ) ~ 1 and the Fourier transform of 
the potential at vanishing wave vector. For a potential of 
size R this leads to a behaviour Ad ~ Vb/m ■ In the 
experiments, however, the repulsive potential is rather 
strong Vb S> /i and, in particular, we want to investigate 
the hard core limit Vb — *■ oo, where perturbation theory 
definitely fails. In this regime the obstacle is equivalent 
to the boundary condition w-n|_R = of vanishing normal 
velocity at \x\ = R. Provided the backflow has the form 
given in (|^) for arbitrary values of r, this immediately 
fixes Ad = R d /(d— 1). With this approximation, (|J) is 
just the flow pattern of an incompressible, ideal classical 
fluid around an obstacle, which has a vanishing drag force 
in a stationary situation v = const. ||. To account for 
the depletion in density due to the finite compressibility 
of a real BEC, we use the following Ansatz 



f(r) = 



— tanh ( r=« 
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for r > R 
for r < R , 
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with £ = h/y/2mfj,, = w2m(VB — A*)/^ 2 an d obstacle 
potential Vb- Equation (m) combines the standard one 
dimensional solution of the GPE at an infinite potential 
step with a one-particle wave function which decays ex- 
ponentially below the barrier Vb ■ The latter is a solution 
of the GPE without the <I> 3 term and thus is valid in the 
limit Vb S> where the external potential due to the 
obstacle is much larger than the mean field interaction 
energy. Matching the logarithmic derivative of f(r) at 
r = R, leads to the following expressions for 5 and A 



A 
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Equation (0) smoothly interpolates on the scale £ be- 
tween the bulk condensate density n = fj,/g far from 
the obstacle to essentially zero at its center r — (note 
R 3> £)• The Ansatz (Q) thus describes a "soft" , penetra- 
ble obstacle as realized experimentally. Although ([?]) is 
not an exact solution of the 3-dimensional GPE, compar- 
ison to numerical solutions shows that the approximation 
is excellent [p^| . 

To calculate the critical velocity of the superfluid flow, 
we follow a method similar to that used to determine the 
stability of electrical supercurrents in a wire by Langer 
and Ambegaokar [jl6| . We consider small deviations 
around the approximate stationary point $ = fe ltp de- 
scribed by equations (g) and (Q) which are conveniently 
written in the form 



$ = (f(r) + w(r)) exp(itp(r)) . 



(10) 



For a general complex function w(r), this Ansatz includes 
both amplitude and phase fluctuations around the sta- 
tionary solution with phase <p(r). Since the potential 
flow around the obstacle described by $ is a stationary 
point of the energy functional ([!]), the expansion of J-[&\ 
up to second order in w 

^[*] + Q^{w} + ... , (11) 

yields a quadratic form Q $ [w] . The eigenvalues of [w] 
give the characteristic curvature of at $ = fe llp in 
function space. In order for T to have a local minimum 
at $, which means that the solution is stable, all the 
eigenvalues of Q<j,[w] must be positive. The lowest ve- 
locity v for which Q has an eigenfunction with negative 
eigenvalue thus gives the critical velocity v c for the onset 
of dissipation in the BEC (l9). The eigenvalue equation 
4%- = Xw can be written in the form 

ow 
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Similar to the situation discussed in |16[ ], it turns out 
that a purely complex w yields the smallest eigenvalue 
and thus the smallest critical velocity. Physically this 
means that the instability is driven by phase fluctuations 
associated with vortex generation, as also found numer- 
ically [0. Using w — iw, the eigenvalue equation ([IT 
can be written as a Schrodinger equation 



2m 
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(13) 



with an effective potential which is determined by in- 
serting the approximations (||) and (^) in equation (|T^). 
The effective potential thus has the following form in d 
dimensions 



Vef fir, i 



Ad 

,.<! 



{(d 2 -2d) cos 2 6» + l} 



^{dcos 2 9-1} 



+ 9f 2 (r). (14) 
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Examining the effective potential, we find that an attrac- 
tive well is formed around 9 — ±7r/2, where bound states 
with negative eigenvalue may develop (see Fig. |l|). The 
depth of this potential well increases with the velocity v, 
and thus there is a critical velocity v c beyond which equa- 
tion ( |l3| ) has a negative eigenvalue. In agreement with 




FIG. 1: Effective potential V eff (r,6 = n/2,v = 0.7c 3 ) for 
Vb = and R = 20£. The solid line is the potential in 
d = 3 and the dashed line in d = 2. 



numerical simulations ]l7j |, 9 = ±ir/2 is the position at 
the side of the obstacle where vortices are generated at 
the breakdown of superfluidity, i.e. it is the barrier for 
vortex creation which vanishes at v c . 

In principle, the velocity at which the lowest eigenvalue 
of the Schrodinger equation with the anisotropic poten- 
tial (|lj) becomes negative may be found numerically. In 
order to obtain an analytical result which displays the 
relevant parameter dependence of the critical velocity, 
however, it is more convenient to use a simple semiclassi- 
cal estimate for the number of eigenstates in the effective 
potential p| 



N sc (v) 
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where Vd is the volume of the unit sphere in d dimen- 
sions. Equation (|l5|) gives the number of states which 
have negative eigenvalue and is just the classical phase 
space volume of the region where V e ff < in units of 
(2irK) d . The condition N sc (v c ) = 1 thus gives a simple 
semiclassical approximation for the critical velocity. 

Introducing reduced units R = R/£ and Vb/p, the 
condition N sc — 1 takes the following form for the cylin- 
drical obstacle 
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For the spherical perturbation the corresponding condi- 
tion for v r reads 
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In (|T^, 17) / is the function ([?]) in reduced units 

ftanh(4 (p- *-§■)) forp>l 
[^exp(E^-l(p-l)J forp<l. 

Equations (|l6|) and (jl^) are the main results of this 
work, providing a simple criterion for the critical velocity 
of a moving cylinder and sphere. Evidently the ratio 
v c /c s depends only on the two dimensionless parameters 
V B /n and Rl£. In the limit R/£ > 1 and V B /p > 1, the 
function ( |l8|) converges to the step function 9(p — 1). In 
this limit the integrals in ( )l6| ) and (|l^) are independent 
of R/£ and Vb / P- The dependence of the critical velocity 
on the obstacle radius is therefore simply 



h 

mR ' 



(19) 



The proportionality factor for the cylinder and sphere 
can be easily evaluated giving v c = 7.6lh/mR and 
v c = A.82h/mR respectively. The precise numerical fac- 
tor is in agreement with the lines shown in Fig. ^[ It 
will change in an exact solution of the Schrodinger equa- 
tion (13), without however, affecting the basic result 
v c = const. h/mR. As shown in Fig. [2] the numerical 
integration of (16[ 17]) gives critical velocities which very 
closely obey the l/i?-scaling ( p~9| ) for radii larger than 
about 10£ or 15£ for the sphere or the cylinder. For 
small radii R < 7£, the critical velocity becomes larger 
than the sound velocity. In this regime our approxima- 
tions are no longer valid and a more microscopic theory 
is necessary to determine v c . Choosing different values 
for Vb/p, it turns out that the result shown in Fig. || 
remains unchanged down to Vb ~ 5/i, beyond which v c 
starts to increase, in agreement with the behaviour found 
numerically || . Numerical work on a moving cylinder in 
a homogeneous BEC fjj observed a critical velocity of 
about 0.45c s for R > £. No dependence of the results on 
the obstacle radius was reported, however. Simulations 
for a trapped BEC and an obstacle with Gaussian object 
potential reported a critical velocity between 0.55c s and 
0.3c s |^|, close to the typical values for v c found in this 
work. 

In the experiments j|, |) the blue detuned laser beam 
which served as a cylindrical obstacle had a radius of 
R w 20£ and a potential barrier Vb — (6.4 — 8)p. The 
measured critical velocity in the trapped BEC was about 
(0.1 — 0.25)c s . For the values of ref. [Q our approach 
yields a critical velocity v c = 0.51c s for a perturbation 
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FIG. 2: Critical velocity v c as a function of the obstacle radius 
for a potential Vb = 10/i. The upper curve presents the 
moving cylinder and the lower one the moving sphere. The 
results obtained from (jlfj) and ([l?]) are shown as crosses, the 
lines exhibit the 1/R dependence. 



moving uniformly in a homogeneous condensate. Com- 
pared with the experiments, this is about a factor of two 



larger. Now, as has been found in recent numerical sim- 
ulations of the GPE H , the measured critical velocity is 
lower than in the bulk for two reasons: First of all, in the 
trapped BEC, the obstacle touches regions of the conden- 
sate with lower density and thus a lower velocity of sound. 
Accordingly vortices will appear first in the outer regions 
of the condensate and penetrate towards the trap center. 
Secondly, due to the oscillatory movement, the obstacle 
creates its own wake, thus lowering the critical velocity 
further. In the experiment || the critical velocity was 
measured for different condensate densities n at fixed R. 
In contrast to our result ([l9j), the critical velocity seems 
to scale with c s ~ \fn or even more strongly with den- 
sity. It is possible that this discrepancy is caused by our 
incompressible fluid approximation for the strength Ad 
of the dipolar backflow. Clearly a more microscopic cal- 
culation of Ad for strong scattering is required to clarify 
the situation and also experiments, in which the effective 
obstacle size R is varied at fixed density. 

In conclusion, we have derived a simple analytical re- 
sult for the critical velocity of a macroscopic moving ob- 
stacle in a BEC. In particular it has been shown that the 
critical velocity scales inversely with the obstacle size, a 
prediction which should easily be checked experimentally. 

We thank W. Ketterle for helpful comments. 
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